Membrane lateral structure: The influence of immobilized particles on domain size 
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In experiments on model membranes, a formation of large domains of different lipid composition 
is readily observed. However, no such phase separation is observed in the membranes of intact 
cells. Instead, a structure of small transient inhomogeneities called lipid rafts are expected in these 
systems. One of the numerous attempts to explain small domains refers to the coupling of the mem- 
brane to its surroundings, which leads to the immobilization of some of the membrane molecules. 
These immobilized molecules then act as static obstacles for the remaining mobile ones. We present 
detailed Molecular Dynamics simulations demonstrating that this can indeed account for small do- 
mains. This confirms previous Monte Carlo studies based on simplified models. Furthermore, by 
directly comparing domain structures obtained using Molecular Dynamics to Monte Carlo simu- 
lations of the Ising model, we demonstrate that domain formation in the presence of obstacles is 
remarkably insensitive to the details of the molecular interactions. 



I. INTRODUCTION 

Membranes are two-dimensional (2D) fluid environ- 
ments, consisting of many different types of lipids and 
proteins [l]. The membrane constituents are not ar- 
ranged randomly, but are expected to be spatially or- 
ganized into domains of different size, composition, and 
dynamics [2]-[6]. In biological membranes this is impor- 
tant because it links to key processes in cells, such as 
signaling, endocytosis, and adhesion [71^. In artificial 
membranes, domain formation is relevant for applica- 
tions, ranging from photolithographic patterning, spa- 
tial addressing, microcontact printing, and microfluidic 
patterning [9l [10]. To identify the factors that control 
domain formation, and to understand their underlying 
physical mechanisms, is therefore of key practical impor- 
tance. 

One challenge is to account for a heterogeneous domain 
structure in thermal equilibrium [2] . Due to line tension, 
one would naively expect a multi-domain structure to 
be unstable as the free energy would be reduced by do- 
main coalescence. Indeed, in model membranes, this is 
precisely what is observed [11] [12]. At high tempera- 
ture, these systems are in a mixed state, but they phase 
separate into (macroscopically large) liquid-ordered and 
liquid-disordered domains at low temperature. At the 
temperature where the phase separation begins to occur 
(the so-called critical temperature), the transition passes 
through a critical point [Til [131 [14] . the critical point, 
the membrane is highly heterogeneous, featuring domains 
of all sizes. Hence, at least in model membranes, a non- 
trivial heterogeneous domain structure arises only near 
the critical point. In line with this, it has been proposed 
that critical behavior might also explain a heterogeneous 
domain structure in biological membranes [15 . We note 
that critical behavior is not the only feasible mechanism. 
Also the coupling between the elastic properties of the 



membrane and its composition has been identified to be a 
key factor affecting domain shape and size [T6H2T1 . In ad- 
dition, the presence of hybrid lipids that collect at the in- 
terface between liquid-ordered and liquid-disordered do- 
mains could also induce a heterogeneous domain struc- 
ture [22]. 

A. Particle immobilization in membranes 

The above mechanisms share in common that they do 
not require the membrane to be coupled to its environ- 
ment in any specific way. However, an essential differ- 
ence between model membranes and real cells may well 
be the fact that the latter are coupled to their surround- 
ings. For example, the cytoskeleton of a cell may induce 
a tension on the membrane affecting its domain struc- 
ture [24 . Another possible effect - the one that we fo- 
cus on in this work - is that the coupling of the mem- 
brane to its environment causes some of its components 
to become immobilized. A possible mechanism leading 
to particle immobilization in a real cell may be an actin 
network underlying the membrane. The actin network 
is attached to the membrane via anchoring proteins and, 
since the actin network is disordered, their positions will 
be random. The anchoring proteins are essentially im- 
mobile compared to the freely diffusing membrane com- 
ponents, and thus act as obstacles for the mobile com- 
ponents [fig.jlja)]. In addition to the anchoring proteins 
themselves, also other molecules can attach to the actin 
strands, which thereby also become obstacles. In fact, the 
distance between obstacles anchored to the actin strands 
can become very small d ~ 2 — 9 nm [23] . 

The immobilization of molecules can also arise in 
vitro in supported membranes [9, 25 -30 . Fluorescence 
measurements have revealed that surface roughness pro- 
foundly limits lipid diffusion, even leading to complete 
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FIG. 1. Schematic sketches (not to scale) showing possible origins of quenched disorder in membranes, (a) In biological 
membranes, quenched disorder may arise from the anchoring of proteins and other molecules to the cytoskeleton network. 
Such anchored proteins are effectively immobilized. The typical distance between immobilized particles can become very small 
d^2 — 9 nm [23]. (b) An analogous in vitro example is provided by supported membranes. In this case, particle immobilization 
may arise from surface roughness or, in the case of larger molecules, from surface friction, (c) In both examples (a) and (b), a 
top- view of the membrane bilayer qualitatively resembles a 2D array of immobilized obstacles (Q), through which the mobile 
particles (M) diffuse. 



immobilization. For example, the immobilized lipid frac- 
tion is around 15% on alumina substrates, and 5% on 
silica Furthermore, due to surface friction, embed- 
ded proteins may also become immobilized [28] [fig.[TJb)]. 
Hence, both in biological membranes as well as in arti- 
ficially supported lipid bilayers, lateral diffusion may be 
envisioned as taking place inside a random network of 
static obstacles [fig. [ijc)]. In physical terms, the pres- 
ence of a random network of static obstacles constitutes 
a form of quenched disorder. 



directly connects to the fundamentals of fluid phase be- 
havior in quenched porous media, for which numerous 
theoretical results are available [44H46] . In particular, 
when conditions (1) and (2) above are met, the obstacles 
induce a change in the universality class, from Ising to- 
ward random-field Ising [47 . As is well known, the latter 
does not support macroscopic domain formation in two 
dimensions [48", "W. Hence, based on the fundamentals 
of statistical physics, particle immobilization provides a 
robust mechanism to account for a heterogeneous mem- 
brane domain structure. 



B. Quenched disorder in membranes: 
Summary of known results 

The above examples show that the presence of 
quenched disorder in membranes is inevitable in many 
practical situations. This naturally raises the question to 
what extent this can account for the properties of a mem- 
brane. Regarding single particle diffusion, it is known 
that lateral diffusion constants are drastically reduced in 
the presence of obstacles 01] [32] . The diffusion constant 
of band 3 in mouse erythrocytes without cytoskeleton, 
i.e. in the absence of quenched disorder, is 50 times larger 
compared to that in healthy cells [32 . Additionally, the 
diffusion becomes anomalous (non-Brownian) on inter- 
mediate time-scales. These results can be explained by 
assuming that the quenched obstacles divide the mem- 
brane into compartments [71 [23l [SU I33H38] . Within a 
single compartment, the diffusion is fast and Brownian, 
while diffusion between compartments is governed by a 
much slower "hopping" dynamics [31 . 

In addition to single particle diffusion, quenched ob- 
stacles also affect the collective behavior of the mobile 
constituents. By collective behavior we mean in this con- 
text the size and shape of the liquid-ordered and liquid- 
disordered domains that can form in the membrane. Pre- 
cisely this question was addressed in a series of recent 
Monte Carlo (MC) simulations [39H43] . A remarkable 
finding is that, provided (1) the quenched obstacles are 
randomly distributed, and (2) show a preferred attrac- 
tion to one of the mobile components, macroscopic do- 
main formation is entirely prevented [41 . This behavior 



C. Aim of paper 

In this paper, we put the recent MC findings [41]- 
^43^ pertaining to particle immobilization and its effect 
on membrane domain formation to a stringent test. It 
should be noted that these previous studies are based on 
rather simple models: all exclude membrane height fluc- 
tuations, most are defined on a lattice, lipid tail degrees of 
freedom are either ignored or only very crudely modeled, 
while solvent and cholesterol are absent. In addition, 
since all deployed the MC method, the dynamics could 
only be approximated or was omitted altogether. Hence, 
it is important to verify if the observed trends in the 
MC studies also survive in more realistic membrane mod- 
els. The aim of this paper is to provide this verification. 
We present a high-resolution Molecular Dynamics (MD) 
study of a phase-separating DPPC/DLiPC/cholesterol 
lipid bilayer in the presence of an immobilized compo- 
nent, which captures important details left out in earlier 
works. 

The outline of our paper is as follows: In Section |ll[ 
we present our MD simulation results, and show that the 
presence of quenched obstacles in a membrane indeed 
results in a heterogeneous structure of small domains. 
In Section |III[ we then compare our MD results to the 
predictions of a simple MC model (the 2D Ising model), 
and demonstrate that the latter is already well-suited to 
describe domain structures. We end with a summary and 
conclusion in Section HVl 



3 



Run 


Type of obstacles 


^max 


I 


None 




25 


/iS 


Ila 


Same on both sides, mobile 


in z-direction 


18 


/iS 


lib 


Same on both sides, mobile 


in ^-direction 


22 


/iS 


III 


On one side only, pinned in 


z-direction 


22 


/iS 



TABLE L Overview of the four MD simulation runs per- 
formed in this work, listing the immobilized obstacle type, 
and the total simulation time tmax- Runs Ila and lib use the 
same obstacle configuration, but different random numbers 
for the initialization of the mobile lipids were used. Note that 
diffusion rates in the MARTINI force field are ^ 4 times larger 
compared to atomistic simulations and experiments [54] . Val- 
ues for time given in this work refer to the simulation time 
coordinate, not to the rescaled (i.e. larger) "physical time" . 



II. MOLECULAR DYNAMICS SIMULATIONS 

A. MD simulation setup 

For the MD simulations presented here, a setup from 
a previous work [50] was used, adapted to our needs by 
changing the size and shape of the membrane patch and 
adding quenched disorder to the system. In the setup 
without disorder, which is simulated first for compari- 
son, the membrane consists of 1532 cholesterol molecules 
and two species of phospholipids, namely 1408 molecules 
of dipalmitoyl-phosphatidylcholine with fully saturated 
carbon tails (DPPC), and 2184 molecules of the polyun- 
saturated phospholipid dilinoleyl-phosphatidylcholine 
(DLiPC). This composition lies deep within the coex- 
istence region of the phase diagram ^5TH53] . Hence, our 
lipids have a strong tendency to demix. 

In the initial state, the lipids form a fiat bilayer patch 
in the (x^y) plane, using a periodic simulation box of 
size 45 X 30 X 11 nm with the same number of DPPC and 
DLiPC lipids in both layers. Within each layer, the phos- 
pholipids are initially randomly mixed. The molecules 
are modeled with the MARTINI force field [54], which pro- 
vides a near-atomic resolution of 3-4 non-hydrogen atoms 
per simulation bead. In order to capture the effect of the 
solvent, 74822 solvent beads were included. The simula- 
tion has been performed with GROMACS [55^ ,56^ Version 4 
at temperature 295 K and a pressure of 1 bar (for further 
details we refer the reader to the original setup [50 ). In 
total, four different simulation runs have been performed, 
summarized in Table [ij and to be discussed further in the 
following sections. 

B. The free membrane 

We first consider the "free" membrane (run I in Ta- 
ble |l|, by which we mean a membrane without any 
quenched disorder. In the absence of quenched disorder, 
model ternary membranes of phospholipids and choles- 
terol readily phase separate into macroscopic (> /im) 




FIG. 2. Membrane without quenched disorder (run I) at 
t = 5 /is, showing the upper leaflet. Due to inter-leaflet cou- 
pling, the domain structure on the lower leaflet is almost iden- 
tical. From the initially mixed state, two domains of DPPC 
(dark blue) and cholesterol (yellow) have formed in a back- 
ground of DLiPC (bright brown). The DPPC domain in the 
foreground has assumed the stripe shape, which is the simu- 
lation analogue of large round domains seen in experiments. 



domains [12]. The domains observed in experiments are 
typically round, as this shape minimizes the length of the 
line interface. The key difference with computer simula- 
tions using periodic boundary conditions is that the equi- 
librium domain shape is not necessarily round. This is a 
crucial point in our analysis and to avoid possible confu- 
sion with related studies [43l |57j we explain it explicitly. 
In a 2D system with periodic boundary conditions, the 
two candidate domain shapes are (1) a circle, or (2) a 
stripe that spans the system along the shorter edge of the 
simulation box [58l [59] . The shape that prevails is the 
one with the shortest line interface, which depends on the 
area fraction of the two coexisting phases and the aspect 
ratio of the simulation box. For the simulations presented 
here, the hallmark of macroscopic domain formation in 
a finite simulation box with periodic boundaries is the 
appearance of the latter stripe geometry. Consequently, 
in our MD simulations, the appearance of such a stripe 
domain will be regarded as evidence for macroscopic do- 
main formation. 

We started the MD simulation of the free membrane 
with a mixed configuration. During the first few /is of 
simulating the membrane, we observe the characteris- 
tics of binary demixing, namely the formation, coarsen- 
ing, and merging of domains. These domains reflect the 
liquid-ordered phase (rich in DPPC and cholesterol) and 
the liquid-disordered phase (rich in DLiPC). Due to inter- 
leafiet coupling [50 , domains in both leaflets arrange in 
similar shapes and at the same locations. After 5 /is, one 
of the domains has assumed the stripe geometry, which 
remains stable for the rest of the simulation [fig. [2]. Fol- 
lowing the explanation above, this indicates macroscopic 
demixing into large DPPC-rich and DLiPC-rich domains, 
such as seen experimentally in free-standing giant unil- 
amellar vesicles [12] [60]. Ultimately, the small DPPC 
domain on the right of fig. [2] will merge with the stripe. 
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FIG. 3. "Early time" domain structures after t = 0.5 /is 
(a), t = 2 /is (b), and t = 5 /is (c) as observed in simulation 
run Ila (cf. Table The obstacles are shown as white; the 
remaining color-coding is the same as in fig. [2] 



but this process is very slow and beyond our computa- 
tional resources. 



C. Membranes with quenched disorder 

To demonstrate how quenched disorder affects struc- 
ture formation in membranes, obstacles have been added 
to the free membrane of Section II B[ Contrary to previ- 
ous studies of membranes with quenched disorder j39ti43] 
that only took into account a flat monolayer our model 
takes into account the full bilayer and allows for mem- 
brane undulations. Hence, there are two additional de- 
grees of freedom we can investigate. First of all, the 
obstacles may be allowed to follow the membrane undu- 
lations in the z-direction, or be fixed in this direction. 
Secondly, the obstacles may be placed at the same lat- 
eral locations in both leaflets, which may be conceived 
to describe a trans-membrane protein, or they could be 
chosen to reside on only one leaflet. 



1. Obstacles that follow membrane fluctuations 

We first consider the case whereby the obstacles have 
been placed at the same locations in both leaflets and 
are allowed to move freely in the 2:-direction. In this sce- 
nario, the obstacles correspond to trans-membrane pro- 
teins that are free to move in the z-direction themselves, 
or so long that the membrane can fluctuate around them. 
To model the obstacles, we chose to take DPPC molecules 
that are not allowed to diffuse laterally. Such obstacles, 
which we refer to as FPPC molecules in what follows, 
naturally attract the DPPC-rich phase. The motivation 
to use DPPC obstacles stems from experimental consid- 
erations: cytoskeleton anchor proteins typically favor the 
liquid-ordered phase [43l[6T|. In principle, we could have 
chosen DLiPC-affine obstacles or a mix of DPPC-affine 
and DLiPC-affine obstacles. Similar results as those pre- 
sented here would be expected in these cases also. 

To generate an initial state of a membrane with obsta- 



cles, we take the initial state from Section II B and chose 



79 random lateral locations in the membrane. Then, in 
each leaflet the 79 nearest phosphate atoms belonging to 
a DPPC molecule are identified and moved to the target 
positions by a series of small translations and subsequent 
energy minimizations (so that the DPPC molecules re- 
main intact and that there is no overlap between lipids). 
Finally, the so-moved DPPC molecules are labeled FPPC 
and the (x, y)-coordinates of their phosphate atoms are 
bound to their target locations via a hard harmonic po- 
tential. During the course of the simulation, the lat- 
eral positions to which the phosphates of the FPPC are 
bound scale with the volume changes of the simulation 
box. Otherwise they are static, effectively making the 
FPPC quenched. 

Two simulation runs with the same FPPC obstacle 
configuration have been performed (Ila and lib of Ta- 
ble|l|, both starting from a state of "mixed" mobile lipids. 
In both simulations, a formation and coarsening of small 
DPPC-rich domains is seen during the first few microsec- 
onds (as shown in fig. [s] for run Ila). After approximately 
5 /is, however, domain growth and merging of domains 
into larger ones is arrested. Beyond this time, the do- 
mains exhibit shape fluctuations, but do not significantly 
change in size or location for the rest of the simulation. 
Most importantly, as shown in fig. [4j neither of the simu- 
lation runs reveals the stripe geometry seen in fig. [2] that 
would indicate domains larger than the size of the simula- 
tion box. We conclude that the obstacles indeed prevent 
macroscopic demixing, and replace it with the formation 
of small domains that appear stable over times at least 
in the order of ^ 10 /is. 

In recent MC studies, it was shown that quenched 
obstacles induce regions in the membrane where ei- 
ther a liquid-ordered or a liquid-disordered domain is 
preferred[41, 42 . At zero temperature this would es- 
sentially fix the domain structure into the groundstate 
configuration, while at finite temperature thermal fluctu- 
ations still allow some freedom as to where the domains 



FIG. 4. "Late time" domain structures of runs Ila (a) and lib (b) at t = 10 /iS. FPPC molecules are shown as white and are at 
the same locations in both systems (the slightly different positions of the white beads between the snapshots result from only 
the phosphate atom of the FPPC being quenched). Otherwise, the color-coding is the same as in fig. [2] In both simulations, 
the stripe domain characterizing macroscopic domain formation is absent. Instead, structures consisting of small domains are 
seen. 




the membrane height fluctuations (run III). Such obsta- 
cles thus locally "pin" the membrane height, which could 
mimic a cytoskeleton anchoring site or, in the case of an 
adhered membrane, a (stiff) ligand-receptor bond [62]. 
This scenario is physically interesting because it probes 
to what extent inter-leaflet coupling transmits the influ- 
ence of quenched disorder from one leaflet to the other. 
To address this scenario in our MD simulations, we per- 



form the same setup steps as in Section II C 1 with the 
same 79 target locations for the obstacles, but we place 
the obstacles on only one of the leaflets and also bind the 
z-coordinate of the FPPC phosphate atoms to the hard 
harmonic potential. 



FIG. 5. Top-view of a membrane with FPPC obstacles in 
only one bilayer leaflet (run III) at t = 10 /is. For clarity, 
this figure does not show the DLiPC or cholesterol molecules. 
DPPC belonging to the leafiet with the FPPC are drawn in 
blue, those of the "free" leaflet in yellow. The domains form 
"in registry" in both leaflets, but they are larger compared to 
the case where FPPC resides in both leaflets. The structure 
is dominated by a single large domain in both leaflets, but it 
is distinctively different from the stripe geometry that would 
indicate macroscopic domains. 



form. This explains why, despite identical FPPC loca- 
tions, the domain structures in fig. |4]^a) and (b) are not 



the same. We return to this point in detail in Section III 



2. Obstacles on a single leaflet 

We now consider the case where the FPPC obstacles 
reside in only one of the bilayer leafiets, and do not follow 



Our simulations indicate that inter-leafiet coupling is 
indeed an important mechanism in domain formation, 
which is consistent with previous studies [171 ISS] • Dur- 
ing the first few microseconds, the formation, coarsen- 
ing, and merging of domains is observed in both leafiets. 
The domains form "in registry" in both leafiets, occur- 
ring at the same (lateral) positions. After 10 /is, most 
of the DPPC molecules are contained in a single domain 
as shown in fig. [5] We thus obtain larger domains com- 
pared to the case of FPPC residing in both leafiets. At 
the same time, we do not obtain the stripe geometry, sug- 
gesting that the domains are not macroscopic. Fig.[5]thus 
presents a "hybrid" structure, between the stripe domain 
of fig. [2] and the small domains of fig. [4] This result shows 
that the effect of quenched disorder is transmitted from 
the leaflet with FPPC to the one without, but also indi- 
cates that the tendency of the leaflet without FPPC to 
form large domains "back-transmits". The net result is 
a structure of flnite-sized domains (i.e. not macroscopic) 
but of a typical size that exceeds the case where FPPC 
resides in both leaflets. 
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FIG. 6. Ratio of the interface length / and domain area 
A of the largest domain for the free membrane (run I), the 
membrane with FPPC in both leaflets (runs Ila and lib), and 
for a membrane with FPPC in only one leaflet (run III). For 
each run, results are shown for each leaflet separately (hence 
8 curves in total). Already at early times, I /A for the free 
membrane differs from the membranes with quenched disor- 
der, and saturates at a significantly lower value. All curves 
shown are interpolations of the raw simulation data. For clar- 
ity, the raw data are only shown explicitly for one curve. The 
scatter in the raw data of the other curves is comparable. 



D. Interface properties of domains 

We now consider a more quantitative measure to dis- 
tinguish membranes with quenched disorder from those 
without. To this end, we consider the ratio I/A of the 
domain interface length / to domain area A. In the 
free membrane, the formation of macroscopic domains is 
driven by the desire of the system to minimize the inter- 
face length, yielding a low value of I /A. In the membrane 
with quenched disorder, where stable small domains are 
expected, I /A will be larger. To measure //A, a 153 x 102 
grid of rectangles is superimposed on the membrane. The 
positions of the phosphate atoms of the DPPC molecules 
are then projected onto the grid. If a grid cell contains 
the phosphate atom of a DPPC molecule, then that cell 
is considered to be part of a DPPC-rich domain. Cells 
with less than four neighbors belonging to a DPPC-rich 
domain are defined as interface. 

In fig. [6j the ratio of the number of interface squares / 
to the number of total squares A for the largest DPPC- 
rich domains are shown for both leaflets and for all our 
simulations. In all cases, the curves first show a rapid 
drop resulting from the initial formation of domains, fol- 
lowed by a saturation as the system equilibrates. The 
saturation values for the free membrane, however, are 
distinctly below those of the membranes with quenched 
disorder. This is consistent with macroscopic domain 
formation in the former, and stable finite domains in the 
latter. Note also that, for the membrane with FPPC ob- 
stacles in only one leaflet, I/A is largest in the FPPC 



containing leaflet. 

E. Choice of order parameter 

In our analysis, domains were identified by their phos- 
pholipid content. An alternative choice is to identify do- 
mains via the tail ordering of the phospholipids, i.e. the 
identification of liquid-ordered and liquid-disordered do- 
mains [57 . To this end, we have also quantified the elon- 
gation of the phospholipid tails via a nematic order pa- 
rameter. DPPC and DLiPC lipids typically have high 
and low nematic order parameters, respectively. Simu- 
lation snapshots in which the phospholipids are color- 
coded according to their nematic order parameter reveal 
the same domain structures as those presented here. The 
identification of domains according to lipid content or tail 
conformation is therefore equivalent, and will not be ex- 
plicitly presented here. 

III. COMPARISON TO THE 2D ISING MODEL 

In Section |TI| we have shown using a detailed mem- 
brane model that quenched obstacles prevent the forma- 
tion of large (macroscopic) domains. This confirms the 
validity of the trends observed in earlier MC studies on 
simple membrane models [4QH43] . One issue that remains 
open is whether the apparent elimination of macroscopic 
domains really is the correct equilibrium behavior, or 
merely caused by the obstacles slowing down the dynam- 
ics of domain formation such that the relevant timescales 
become inaccessible in the simulation. Because of their 
high computational demancQour MD simulations alone 
cannot resolve this issue. However, in MC simulations of 
simple models, where the proper equilibrium behavior (of 
these simple models) can be accessec|^ it is predicted that 
macroscopic domain formation is not merely delayed, but 
in fact eliminated (the theoretical motivation is the Imry- 
Ma argument [41, 48 ). In this section, we will show that 
these MC models make further predictions that can be 
tested against our MD results. Provided one uses the 
same positions for the obstacles, the resulting domain 
structure is remarkably insensitive to the model details. 
In fact, the minimal model requirement is that the unlike 
lipid species repel each other, which is already contained 
in the 2D Ising model. In what follows, we will show 
that the domain structures obtained in fig. |4] with MD, 
can indeed be reproduced in MC simulations of the 2D 
Ising model. 

In the 2D Ising model, each lattice site i of a 2D square 
lattice is occupied by a spin variable Si = ±1. In the con- 



^ Each of the MD runs of Table ^ required several months of sim- 
ulation with 32 processors in parallel. 

^ In the case of MC simulations of the 2D Ising model, the equiva- 
lent of one MD run takes only a few seconds on a single processor. 
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FIG. 7. Typical domain structures of the "Ising" membrane 
MC simulation at J/k^T — 0.6. (a) Domain structure in the 
absence of quenched obstacles, and with 681 DPPC (blue) 
and 1053 DLiPC lipids (orange). The stripe geometry is re- 
covered, in agreement with the discussion of Section jll B| The 
snapshots (b,c,d) show typical structures in the presence of 78 
FPPC obstacles (red). In each of these snapshots, the FPPC 
configuration is the same, and chosen to reflect the MD con- 
figuration of fig. |4] (runs Ila and lib). The differences in the 
snapshots (b,c,d) reflect the thermal fluctuations. 



text of membranes, Si = +1 (—1) means that site i con- 
tains a DPPC (DLiPC) lipid, while cholesterol is ignored 
on this level [42 . The energy of an Ising spin config- 
uration is given hy E = ~J^<ci j> ^i^j^ with coupling 
constant J > (the sum runs over nearest neighbors). 
We use the Ising model to mimic the domain structures 
of fig. |4] by simulating a 51 x 34 lattice. This approx- 
imates the aspect ratio of the MD simulation box, as 
well as the total number of phospholipids in a single 
leaflet. To represent the FPPC obstacles, the lattice is 
superimposed on the initial state of the MD simulations 
with obstacles. Those sites containing an FPPC phos- 
phate atom subsequently have their spin value frozen to 
Si = +1. The remaining sites are randomly initialized 
with values Si = ±1, under the constraint that the re- 
sulting (DPPC:DLiPC) composition matches that of the 
MD simulation as well as possible. 

We then use MC simulations with Kawasaki dynam- 
ics [6T to determine equilibrium domain structures. 
That is, in each MC move, two non-FPPC lattice sites 
are chosen at random. Then, the corresponding spin 
values are exchanged with the Metropolis probability. 
Pace = min[l,exp(-A£^/A:BT)], with AE the energy dif- 
ference caused by the swap, /cb the Boltzmann constant, 
and T the temperature. In the absence of quenched dis- 
order, the Ising model exhibits macroscopic phase sepa- 
ration when J//cbT > ^ ln(l + V^) ^ 0.44, as was shown 
by Onsager ^65j. Note that J/k^T is the only free pa- 
rameter in the context of this work. 

At J/k^T = 0.6, the parameter we have chosen to 
use for the comparison with the MD results, the system 
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FIG. 8. The ensemble averaged snapshot (i.e. the Ai^s) of the 
"Ising" membrane simulation corresponding to the obstacle 
configuration of fig.[7|^b,c,d) (which, in turn, was derived from 
the MD obstacle configuration of fig.[4|. Preferential locations 
for the appearance of domains are seen, which coincide well 
with the domain structures obtained in the MD simulations 

of fig.m 



therefore assumes the stripe geometry [fig. [7|^a)]. In the 
presence of quenched disorder, macroscopic phase sepa- 
ration is eliminated [fig.[7|^b,c,d)]. Instead, finite domains 
that qualitatively resemble those of fig. [4] are seen. This 
similarity is quite remarkable, given that the simulation 
models used are entirely different (only the positions of 
the obstacles are the same). It shows that the generic 
features of the domain structure are essentially encoded 
in the obstacle configuration, and rather insensitive to 
the details of the interactions. To put this statement on 
a more solid basis, we now consider a quantitative mea- 
sure of the correspondence between the domain struc- 
tures seen in the two models. 



A. Correlation between the models 

Because of the Ising model's computational simplic- 
ity, many different equilibrium states can be generated 
quickly by means of MC simulation. In fact, for the lat- 
tice size used here, a representative set of the full equi- 
librium ensemble can be generated. This allows us to 
create an ensemble averaged snapshot by generating a 
lot of snapshots such as shown in fig. [7|^b,c,d) and evalu- 
ating the average spin value Ai at each site i. As shown 
in fig. [8) this average snapshot reveals structure, which 
is due to the fact that the obstacles induce regions that 
locally prefer one of the lipid species ^2 . With this 
average structure, we demonstrate the correspondence 
between the MD simulation results and those of the Ising 
model in two steps: First, we define a compatibility mea- 
sure C that quantifies how much a given snapshot of ei- 
ther model resembles the average spin values Ai. Then, 
we evaluate the range of compatibilities within the Ising 
model itself, and compatibilities obtained from a struc- 
ture of domains at random locations (i.e. we want to 
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rule out that an observed correspondence is due to pure 
chance). It wih be shown that the compatibihties of our 
MD results fall into the region of the former, and distinc- 
tively disagree with the predictions for the latter. 

We define the compatibility of an Ising model snapshot 
with the Ai as 

i 

where the summation runs over the N lattice sites that do 
not contain an obstacle. This measure quantifies the de- 
gree of correlation between a single snapshot (character- 
ized by spin values Si) and the ensemble-averaged snap- 
shot of fig. [s] (characterized by averages Ai). To calculate 
the compatibility of an MD snapshot with the averaged 
Ising structure of fig.jsj we formally use eq. ([T]), but with 
proper MD equivalents of the appearing terms. We con- 
sider each leaflet of the membrane separately, so the sum 
in eq. (fTl) runs over all DPPC and DLiPC molecules in 
the leaflet, where N is their total number. If lipid i is 
a DPPC (DLiPC) lipid, then the equivalent spin value 
Si = {si = —1). The corresponding value of Ai is de- 
termined by projecting the leaflet onto the (x, y) plane, 
and to then superimpose the averaged Ising structure of 
flg. [8| For each lipid i, the corresponding Ai is taken 
from the grid cell in the averaged structure that contains 
the lipid's center of mass. 

Since the Ai will in general have values between —1 and 
+1, and since typical snapshots [cf. flg. [7|^b,c,d)] do not 
exactly have the domain structure of the average snap- 
shot [flg. [s], we do not expect the compatibilities of our 
MD simulation to reach C ^ 1. Instead, the proper C 
values to compare against are those of the Ising model 
snapshots themselves, which are straightforward to ob- 
tain: We simply let the MC simulation that was used to 
produce flg. |8] run a second time and collect a histogram 
of C values. The histogram is Gaussian, with mean Cav 
and width cr; the upper shaded region in flg. [9] marks the 
interval (Cav ± cr)/C'av In case the domains obtained in 
the MD simulations are compatible with the averaged 
Ising structure of flg. [Sj we expect the corresponding 
compatibilities to be inside or at least in the proxim- 
ity of the upper shaded area. Indeed, as flg. |9] shows, 
the MD compatibilities profoundly increase with time, 
indicating that the corresponding domain structures in- 
creasingly correlate with the equilibrium domains of the 
Ising simulation. As may be expected, the correlation is 
most pronounced in the case where the FPPC obstacles 
reside in both leaflets, where an excellent agreement with 
the Ising model is found. However, also when FPPC re- 
sides in only one leaflet, the correlation remains. In this 
case, the correlation is most pronounced in the leaflet 
containing the obstacles, as might be expected. 

To demonstrate that the agreement between the do- 
main structures of the Ising model and the MD simu- 
lations does not merely occur by chance, we also mea- 
sure the typical compatibilities that one obtains in case 
the domains appear at random locations. To generate 
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FIG. 9. Compatibility of our MD simulations with the aver- 
age domain structure of the Ising model in fig. [S] The com- 
patibility is normalized by Cav, which is the average compat- 
ibility of the Ising model with its own average structure. The 
upper shaded region corresponds to one standard deviation 
around this average. The inner and outer lower shaded areas 
represent one standard deviation regions that are expected 
for randomly located domains (corresponding to r = 2 nm 
and r = 7.5 nm, respectively; see details in text). The MD 
simulations with FPPC obstacles in both leaflets (runs Ila 
and lib) equilibrate to states consistent with the Ising model 
(i.e. around the upper shaded region) and are clearly incom- 
patible with random domains. In the case of FPPC in only 
one leaflet (run III), the C- values are still distinctly differ- 
ent from random domains (i.e. outside the lower shaded re- 
gion), but neither the FPPC leaflet (solid curve) nor the "free 
leaflet" (dashed curve) appear to be fully compatible with the 
Ising model. 



a random domain structure, we take the initial state of 
the MD simulation, but with all phospholipids (except 
FPPC) being DLiPC. Then, we pick a random location 
in the membrane and change the identity of all DLiPC 
lipids within a radius r of this location to DPPC. This 
step is repeated until the (DPPCiDLiPC) ratio equals 
that of the simulation (note that the resulting domains 
are usually not circular, since the circular regions picked 
at each step can overlap). Since this process is fast, we 
can generate many such random domain structures and 
collect a histogram of compatibilities. This histogram is 
also Gaussian, but centered around a lower average value, 
and with a different width that depends on r. The range 
of compatibilities expected for these domain structures 
is marked in fig. [9] by the lower shaded areas. As can be 
seen in the figure, they are clearly distinct from the MD 
curves (which holds for all values of r) . We therefore con- 
clude that the correspondence between domains obtained 
in our MD simulations and the Ising model is not coinci- 
dental. Despite the very different interaction potentials 
of the two models, both yield equivalent domain struc- 
tures. Consequently, we expect that the domains seen in 
our MD simulations are indeed equilibrated states, as op- 
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posed to kinetically trapped states that would eventually 
macroscopically phase separate. 

We emphasize that our findings are not unique to 
J/k^T = 0.6 considered here. We have also performed 
the analysis for J/k^T = 0.45, 0.5, . . . , 0.65, and always 
found the same scenario: A reasonable agreement be- 
tween the MD simulations with obstacles on both sides 
and the Ising model, while the "pure chance" regime is al- 
ways excluded. The reason for this agreement is that over 
a wide range of temperature, the average domain struc- 
ture [fig. [s] qualitatively remains the same, and merely 
becomes sharper as the temperature is decreased . 

IV. SUMMARY 

In this work, we have presented results from MD sim- 
ulations of a DPPC/DLiPC/cholesterol membrane at 
near-atomic resolution in the absence and in the presence 
of quenched disorder. Our results show that static com- 
ponents in the membrane destroy the formation of large 
domains. Instead, a heterogeneous structure of small do- 
mains is seen [fig. [4]. For this to happen, it suffices that 
the obstacles are present in only one of the leafiets, since 
quenched disorder effects can be transmitted to the other 
leaffet via inter-leaffet coupling [fig. [s]. As the results 
have been obtained for a system with a strong demixing 
tendency, the elimination of macroscopic domain forma- 
tion should occur for other lipid compositions also. A 
coupling of a cell to its surroundings or its cytoskeleton 
that causes an immobilization of membrane components 
may therefore explain the appearance of raft-like domains 
in thermal equilibrium. 

By means of a suitably constructed compatibility mea- 



sure we have also shown that when the obstacles reside in 
both leafiets, the domain structure can be predicted with 
very simple models, such as the Ising model (cf. fig. [9|. 
We do not expect this to be a feature unique to the 
Ising model, but instead consider this result a signal that 
structure formation in the presence of quenched disorder 
is universal and does not strongly depend on the mi- 
croscopic details of the system. Other similarly simple 
membrane models [4TJ [57l [66] are expected to be equally 
suited to predict domain structure. The exception is in 
cases where the obstacle configurations differ between the 
leaffets. In this situation, a simple monolayer model may 
be inadequate due to the neglect of inter-leaffet coupling. 
The construction and use of a related bilayer model is 
conceivable but has not been tested. 

If one accepts the evidence that a MC simulation of a 
simple model can adequately describe the equilibrium do- 
main structure of membranes in the presence of quenched 
disorder, their comparably low computational demands 
allow considering questions that are inaccessible in de- 
tailed MD simulations. The possibility to probe larger 
length scales permits the investigation of obstacle con- 
ffgurations with a structure on their own [42l [32, and 
connect to the length scales of experimental visualiza- 
tion methods [671 HI] ? while the option to simulate many 
different setups allows to compare different positionings 
and different types of static obstacles [39^ .41 j. 
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